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ABSTRACT 

A  computer  simulation  model  is  developed  treating  wave 
propagation  in  a  random,  inhomogeneous  ocean  as  transmission 
through  a  linear  time-invariant,  space-variant  random  communi- 
cation channel.   The  ocean  volume  is  modelled  by  an  index  of 
refraction  which  is  decomposed  into  a  depth-dependent  deter- 
ministic part  and  a  depth-independent  Gaussian  zero-mean 
random  part.   Computer  simulated  output  electrical  signals 
were  generated  that  depend  on  the  complex  frequency  spectrum 
of  the  transmitted  electrical  signal,  the  far-field  beam 
pattern  of  the  transmit  array  and  the  random  transfer  function 
of  the  ocean  medium.   Output  was  generated  for  different  test 
cases.   In  all  cases  the  transmit  electrical  signal  was  repre- 
sented by  a  finite  Fourier  series  and  random  cases  were 
modelled  by  a  random  number  generator.   The  computer  simulated 
output  electrical  signals  were  then  processed  by  a  3-D  DFT 
beamformer  and  the  results  for  the  deterministic  inhomogeneous 
and  random  inhomogeneous  cases  were  compared  to  the  homogeneous 
non-random  case  in  order  to  study  the  effects  of  the  medium 
on  signal  distortion  and  source  localization. 
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I.   INTRODUCTION 

Based  on  the  linearized  acoustic  wave  equation  for  small 
amplitude  acoustic  signals,  sound  propagation  through  the  ocean 
medium  can  be  modelled  as  a  signal  passing  through  a  linear 
time-variant,  space-variant  random  filter.   The  term  "space- 
variant"  implies  that  the  sound  speed  profile  is  a  function 
of  position.   The  space  variant  property  results  in  scatter 
or  angular  spread  due  to  refraction.   If  the  filter  is 
"space-variant,"  then  an  isospeed  medium  is  implied.   As  a 
result,  there  will  be  no  refraction,  and  hence,  no  scatter  or 
angular  spread  since  the  sound  rays  will  be  travelling  in 
straight  lines.   The  term  "time-variant"  implies  motion  of, 
for  example,  discrete  point  scatterers,  the  ocean  surface, 
and  the  transmit  and  receive  arrays  (apertures) .   The  time- 
variant  property  results  in  both  Doppler  spread  and  spread 
in  time  delay. 

By  using  a  linear  system's  theory  approach,  different 
propagation  paths,  such  as  direct  paths,  bottom  reflected 
paths  and  surface  reflected  paths,  can  be  modelled  as  the 
outputs  from  linear  filters.   Furthermore,  different  trans- 
mit signals  and  transmit  and  receive  far-field  directivity 
functions  can  easily  be  coupled  to  various  transfer  functions 
of  the  random,  inhomogeneous  ocean  medium  in  order  to  study 
their  effects  on  signal  distortion  and  source  localization 
using  various  space-time  signal  processing  algorithms. 


Work  based  on  this  approach  has  been  done  by  Laval  [Refs. 
5,6],  Laval  and  Labasque  [Ref.  7],  Middleton  [Refs.  8,9] 
and  Ziomek  [Ref.  1].   Middleton  [Refs.  8,9]  described  propaga- 
tion phenomena  in  a  random  inhomogeneous  ocean  with  the  use 
of  space-time  operators  but  did  not  concern  himself  directly 
with  the  derivation  of  random,  time-variant,  space-variant 
transfer  functions. 

Ziomek  [Ref.  2]  derived  a  time-invariant  space-variant 
random  transfer  function  of  the  ocean  volume  based  on  the  WKB 
approximation.   Ziomek  [Ref.  3]  also  derived  an  equation  for 
the  random,  output  electrical  signal  at  each  element  in  a 
receive  planar  array  of  complex  weighted  point  sources  in 
terms  of  the  frequency  spectrum  of  the  transmitted  electrical 
signal,  the  transmit  and  receive  arrays,  and  the  random  ocean 
medium  transfer  function. 

The  purpose  of  this  thesis  is  to  implement  on  a  computer 
the  equation  for  the  output  electrical  signal  as  derived  by 
Ziomek  [Ref.  3] .   This  mathematical  computer  model  simulates 
the  effects  of  wave  propagation  through  a  random  inhomogeneous 
ocean.   The  ocean  volume  is  modelled  by  the  transfer  function 
derived  in  Ziomek  [Ref.  2]  and  incorporates  an  index  of 
refraction  which  is  a  function  of  depth.   The  index  of  refrac- 
tion is  decomposed  into  a  depth-dependent  deterministic  com- 
ponent and  a  depth-independent  zero  mean  Gaussian  random 
component . 

Section  II  is  devoted  to  a  theoretical  description  of  the 
computer  model.   The  notation  and  geometry  of  the  computer 
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model  are  explained  and  the  fundamental  equations  on  which 

the  computer  model  is  based,  are  stated.   Also,  some  additional 

constraints  are  made  to  allow  for  simpler  implementation. 

Section  III  describes  the  numerical  techniques  used  and 
the  resulting  mathematical  procedures  which  are  implemented 
in  the  program.   Furthermore,  a  brief  outline  of  the  computer 
program  is  given. 

In  Section  IV  a  number  of  test  cases  are  defined  which 
were  used  to  test  the  computer  model.   Computer  simulated 
output  electrical  signals  were  generated  for  test  cases 
incorporating  deterministic  homogeneous,  deterministic 
inhomogeneous  and  random  inhomogeneous  medium  models.   The 
output  signals  were  processed  by  a  3-D  DFT  beamformer  and 
the  results  for  the  deterministic  homogeneous  case  were  used 
as  the  baseline  case.   Both  deterministic  and  random 
inhomogeneous  cases  were  compared  to  this  baseline  case  in 
order  to  study  the  effects  of  the  medium  on  signal  distortion 
and  source  localization. 
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II.   THEORY  ON  THE  COMPUTER  SIMULATION  MODEL 

A.   MODEL  DESCRIPTION 

The  mathematical  computer  model  calculates  the  output 
electrical  signal  at  each  element  of  a  receive  planar  array 
of  point  sources  in  terms  of  the  complex  frequency  spectrum 
of  the  transmitted  electrical  signal,  the  transmit  and  receive 
arrays,  and  the  random  ocean  medium  transfer  function.   The 
theory  for  the  model  is  based  on  time-variant,  space-variant 
random  filters  as  discussed  in  Ziomek  [Ref.  1] .   Ziomek  [Ref. 
2]  derived  a  time-invariant  space-variant  transfer  function 
for  a  random  ocean  volume  where  the  index  of  refraction,  which 
is  composed  of  a  deterministic  and  a  random  component,  is  a 
function  of  depth.   This  transfer  function  is  based  on  the 
WKB  approximation,  which  is  a  solution  to  the  wave  equation 
when  the  index  of  refraction  is  a  function  of  depth.    Further- 
more, Ziomek  [Ref.  3]  derived  an  equation  for  the  random 
output  electrical  signals  appearing  at  each  element  in  a 
receive  planar  array  of  complex  weighted  point  sources  in 
terms  of  the  frequency  spectrum  of  the  transmitted  electrical 
signal,  the  transmit  and  receive  arrays,  and  the  random  ocean 
transfer  function.   The  computer  model  combines  and  implements 
the  equations  derived  by  Ziomek  [Refs.  2,3]  for  the  transfer 
function  and  the  output  electrical  signals. 

1 .   Geometry  of  the  Model 

Figure  1  shows  the  geometry  of  the  communication 
channel  as  simulated  by  the  computer  model.   The  model 
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describes  the  wave  propagation  between  a  transmit  planar  array 
and  a  receive  planar  array  whose  centers  are  located  at 
(x  ,y  ,  z  )  and  (x  ,y  ,  z  ) ,  respectively.   The  x,  y  and  z  axes 
describe  cross  range,  depth  and  range,  respectively.   The 
transmit  array,  which  is  a  parallel  to  the  XY  plane,  consists 
of  M  xN  complex  weighted  point  sources,  where  M  and  N  are  odd. 
The  sources  are  equally  spaced  in  the  X-  and  Y-directions 
with  spacings  of  d   and  d  meters,  respectively.   Note  that 
in  general  d   ^  d  .   The  complex  weights  are  assumed  to  be 
separable  and  given  by 


c      a  exp  { i  0  }  (2.1) 

m      m  ^    J    m 


for    the  weighting    in   the   X-direction    and  by 


dn      =      bnexp{j^n)  (2.2) 


for  the  weighting  in  the  Y-direction.   The  combined  weight 

factor  for  each  transmit  element  is  c  d  ,  where  m  and  n  are 

m  n 

the  indices  describing  the  transmit  element  at  position 
(x  +md  ,y  +nd  ,z  ) .   The  complex  weights  are  used  for  ampli- 
tude shading  and  beamsteering  of  the  transmit  beam  pattern. 
The  receive  planar  array,  which  is  also  parallel  to  the  XY 
plane,  consists  of  M'  xN'  elements,  where  M1  and  N '  are  odd. 
The  sources  are  equally  spaced  in  the  X-  and  Y-directions  with 

spacings  of  d'  and  d1  meters,  respectively,  where  in  general 
x       y 

d'  ^   d ' .   The  oosition  of  a  receiver  element  with  indices 
x     y 

(i,q)  is  given  by  (x  +id',y  +qd ' , z  ). 
-a      =>       -;    r    x^r^yr 
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The  speed  of  sound  in  the  medium  is  assumed  to  be 
only  a  function  of  depth  and  is  given  by  c (y) .   The  reference 
speed  of  sound  c   is  the  speed  of  sound  at  the  center  of  the 
transmit  array,  i.e.,  c   =  c(y  ).   The  index  of  refraction 
n(y)  is  a  function  of  depth  and  is  composed  of  a  determinis- 
tic component  and  a  random  component: 


n(y)   =   CQ/c(y)   =   nD(y)  +  nR(y)  (2.3) 


where  n   is  the  deterministic  component  and  n_  is  the  random 

component  of  the  index  of  refraction.   The  initial  propagation 

vector  will  be  denoted  by  k  ,  k   =  27rf/c   is  the  magnitude 

■'—no        o  ^ 

of  the  initial  propagation  vector  at  depth  y  =  y   (i.e.,  at 

the  center  of  the  transmit  array) .   The  components  of  the 

initial  propagation  vector  along  the  x,  y  and  z  axes  are 

described  in  terms  of  the  direction  cosines  u  ,  v   and  w  : 

o    o      o 


k    =   k  u  (2.4) 

x      o  o 


k    =   k  v  (2.5) 

y       o  o 


k    =   k  w  (2.6) 

z      o  o 


and 


2  2     2 

w=l-u-v  (2.7) 

o  o     o 


Note  that  the  magnitude  of  the  initial  propagation  vector  at 
the  transmit  elements  which  are  located  at  a  depth  nd   meters 

y 
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from  the  center  of  the  transmit  array  is,  in  general,  not 
equal  to  k  and  is  given  by  k  : 


k    =   27rf/c(y  +nd  )  (2.8) 

n         '  2o        y 


One  can  also  use  spatial  frequencies  (in  cycles/meter)  to 
denote  directions  of  propagation: 


f   =s   u  /A;   f    =   v  /A;   f    =  w  /A       (2.9) 
x       o       y       o  z       o 


2 .   Transfer  Function 

The  medium  transfer  function  H   is  given  by  a  function  of 
frequency  f,  spatial  frequency  f   and  receiver  depth  y  for 
a  given  source  depth  y   [Ref .  2] : 


HM(f'VY)   =   AM(f'fy;^)eXp{j[eMD(f'fy^) 

+  eMR(f/fy;y)]>  (2-io) 


where  A   is  the  amplitude  term  given  by 


am  ■   <2»y~1/2  =   <Vo,_1/2  (2-11) 


and  B.,~,    6AA-  are  the  deterministic  and  random  phase  terms, 
MD    MR  ^ 

respectively : 


2  rY         2 

a    /4-rrf 
MD  o  y 


=       (-k^/47Tf    )  /       [n*U)    -l]dc  (2.12 


^o 
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.    =   (-ko/2lTfy)    /Y  %(5)nRC;)d5  (2.13 


This  transfer  function  is  valid  in  the  absence  of  turning 
points  and  when  the  following  inequality  is  obeyed: 


n2(y)  -  l|/v2   <   1  (2.14 


In  the  case  of  a  speed  of  sound  profile  with  a  constant 
gradient  g,  the  deterministic  component  of  the  index  of 
refraction  becomes: 

nD(y)   =   co/(cq  +g(y-yo))  (2.15) 

Substitution  into  equation  (2.12)  results  in 


k   C 
°-{~[nn(y)  -1]  +  y  -  y  )  (2.16 


MD      2v  L  g  L  D  J '  J    J  J o 

o  y 

As  far  as  the  random  phase  term  is  concerned  we  will 
proceed  by  assuming  that  the  random  component  of  the  index  of 

refraction  is  not  a  function  of  depth  and  can  be  represented 

2 
by  a  Gaussian,  zero  mean  random  variable  with  variance  a  : 


nR(C)   =   nR  =   anNR  (2.17) 

2 
where  n^  is  N(0.a  )  and  n>TT^,  which  is  the  normalized  random 
R  NR 

component  of  the  index  of  refraction,  is  N(0,1) .   When  n 
is  not  a  function  of  depth,  equation  (2.13)  becomes: 


17 


'MR   -   -  ^  "R    lY    nDU)dC  <2-18) 

y0 


Substituting  the  equation  of  n   for  a  constant  gradient  into 
equation  (2.18)  leads  to: 

k   c 

8m   =   —  —  n_  ln{nn(y)  }  (2.19) 

MR      v   g   R      D  d 

o  3 

For  a  homogeneous  medium,  equations  (2.16)  and  (2.19)  become: 

6MD   =   0  (2.20) 

6MR   "   "  ^  V*"^  (2-21» 

o 

Equations    (2.11)     through    (2.13)     for    the    amplitude    and 
phase    terms   of   the   transfer    function    show    that   the  medium  has 
the   effect   of   angle   modulating    the   transmitted    sound    field, 
also   referred   to    as    scattering. 
3 .       Output    Electrical    Signal 

The    output   electrical    signal    at    a    receive    location 

(x   +id' ,y   +qd' ,z    )     is    given   by    [Ref.    3]: 
jl  vc       jl  y       jo 

oo  11 

y(tfX  +id^y  +qd',z  )     =         /     X(f)        /         /     f2   • 
y  -oo  -la 

q 

(N-D/2 


J  (d  /c)BL,(f,f  ;y)exp{-jk  v  AY     }exp{-jk    (1-u  -v  ) v   AZ} 

n=-  (N-D/2     n    n    ™        y  FJnoqn       ^     J  n         oo 


(M-D/2  2      2 

V  c  exp{-jk  u  AX.    }dv  du  exp{  j2Trft}df ;     u  +v     <  1        (2.22) 

/,„iw->mno:Lm        oo^J  oo  — 

m=- (M-D/2 
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where 


aq   =   r|n2(yr+qd')  -1111/2 


AY     =   y   -  y   +  qd *  -  nd 
qn      ■*  r    J  o    ^  y      y 


AX.       x   -  x   +  id'  -  nd 
lm       r     o      x      x 


AZ   =   z   -  z 
r     o 


X(f) :   complex  frequency  spectrum  of  the  transmitted 
electrical  signal. 


c    =   c(y  +ndy) ;  speed  of  sound  at  transmit 
element  with  index  n. 


:    =   2irf/c(y  +  ndv)  ;  magnitude  of  initial  propagation 
vector  at  transmit  element  with  index  n. 


c  ,  d  :   complex  weights  of  the  transmit  array 


H  (f , f  ;y) :   Random  time-invariant  space-variant 
"transfer  function  of  the  ocean  medium- 
Note  that  H^  is  a  function  of  both  indices 
q  and  n,  i.e.,  it  is  a  function  of  the  y 
coordinate  of  the  source  location  yQ  +  nd 
and  receiver  location  y  +  qd '  .   This  can^ 
also  be  expressed  by  using  the  notation 
H(k  ,v  ,n,q)  expressing  H   as  a  function  of 
wave  number  k  ,  direction  cosine  v   and 
indices  n  and  q. 


The  above  equations  show  how  the  medium  transfer 
function  is  used  in  describing  the  acoustic  field  at  a  receiver 
location  as  a  function  of  the  acoustic  field  at  a  transmit 
location.   The  far-field  directivity  functions  of  transmit  and 
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receive  arrays  are  used  to  couple  an  electrical  signal  to 

the  medium  and  vice-versa. 

The  reqion  of  inteqration  over  direction  cosine  v 
^  ^  o 

is  limited  by  a  .   This  expresses  the  fact  that  the  transfer 

q 

function  is  not  valid  for  direction  cosine  values  v   approach- 

o 

ing  0  (close  to  a  turning  point) .   In  most  applications  the 
limits  of  integration  for  both  integrals  can  be  further  limited 
as  shown  in  Section  III. 

B.   CONSTRAINTS  ON  THE  MODEL 

In  this  section  the  inherent  constraints  on  the  model  will 
be  investigated.   Also,  some  additional  constraints  will  be  made 
1 .   Geometrical  Constraints 

These  constraints  concern  the  transfer  function  which 
is  only  valid  under  the  two  following  restrictions: 

a.   The  assumption  made  in  the  derivation  of  the  transfer 
function,  in  order  to  validate  a  binomial  expansion  was  [Ref.  2] 


[n2(y)  -U/V2]  |   <  1  (2.23) 


Thus  before  using  the  computer  model  one  should  evalu- 
ate equation  (2.23)  for  the  chosen  sound  speed  profile  and 
array  locations  to  make  sure  the  inequality  is  obeyed. 

b.   Absence  of  turning  points.   The  occurrence  of  turning 
points  in  a  particular  transmitter/receiver  geometry  depends 
on  the  initial  direction  of  propagation,  described  by  the 
angle  3   between  the  initial  propagation  vector  and  the 
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Y-axis.   Using  Snell's  law  we  find  that  at  the  turning  point 
depth  yT: 


c(yT)   =   co/sin(3Q)  (2.24 


The  relation  between  the  angle  3   and  v   is  given  by 


v   =   cos(B  )  (2.25) 

o  o 


Combining  equations  (2.24)  and  (2.25)  into  an  inequality 
which  has  to  be  obeyed  to  avoid  a  turning  point  results  in 


c 


n(y)      ^y   >  A    -    v2q  (2.26) 


and  we  conclude  that  obeying  the  inequality  (2.23)  also 
guarantees  the  absence  of  turning  points. 
2 .   Sound  Speed  Profile 

The  results  presented  in  this  thesis  are  based  upon 
a  sound  speed  profile  which  is  assumed  to  have  a  constant 
gradient  g  and  a  random  component  of  the  index  of  refraction 
nR  which  is  not  a  function  of  depth.   Note  that  the  equations 
derived  in  Ziomek  [Refs.  2,3]  are  more  general  and  that  the 
computer  program  described  in  this  thesis  can  be  easily  modi- 
fied to  handle  a  more  complex  sound  speed  profile.   For  purpose 
of  computer  simulation,  the  random  component  of  the  index  of 

refraction  was  assumed  to  be  a  zero  mean,  Gaussian  random 

2 
number  with  variance  a   and  was  modelled  by  a  random  number 
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generator.   These  restrictions  make  it  possible  to  evaluate 
the  transfer  function  with  the  closed  form  expressions  given 
in  equations  (2.11),  (2.16)  and  (2.19)  for  amplitude,  deter- 
ministic phase  component  and  random  phase  component  of  the 
transfer  function,  respectively.   In  the  calculation  of  the 
random  phase  terms  of  the  transfer  function,  a  different 
random  number  was  drawn  for  each  combination  of  transmit 
element  depth  y  +  nd   and  receive  element  depth  y  + qd ' , 
thereby  assuming  that  the  random  phase  terms  for  the  different 
transmission  paths  are  uncorrelated . 

3 .   Frequency  Spectrum  of  the  Transmitted  Electrical 
Signal 

The  transmitted  electrical  signal  was  represented  by 

a  finite  Fourier  series  with  fundamental  period  T   and  maximum 

r       o 

frequency  f     Hz: 
^     J   max 

f      =   KMAX  •  =i-  (2.27) 

max  T 

o 

where  KMAX  denotes  the  total  number  of  harmonics  in  the  signal, 
This  representation  does  not  impose  a  severe  restriction  on 
the  input  signal,  since  every  finite  energy  signal  can  be 
represented  in  the  sense  of  a  least  mean-square  error  by  a 
finite  Fourier  series.   The  expression  for  the  frequency 
spectrum  becomes : 


KMAX 
X(f)   =       [c,6(f-kf  )  +  c*5(f  +kf  )]      (2.28) 

i -i     K  O        K  O 
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where  f   is  the  fundamental  frequency  1/T   and  the  DC  term 
o  ^     J  o 

is  assumed  to  be  zero.   The  complex  Fourier  series  coefficients 
c,  determine  the  amplitude  and  phase  relations  among  the 
different  frequency  components.   The  integral  in  equation 
(2.22)  now  reduces  to  a  summation  and  the  expression  for  the 
real  output  electrical  signal  at  a  receiver  element  with 
indices  i  and  q  becomes : 


KMAX  ?  1  1 

y(t,i,q)     =     2Re  {   I       c,(kf Y       j         j     • 

k=l  -1       a 

q 

(N_1)/2  2  2     2  1/2 

y  (d  /c  )R,(k   ,v  ,n,q)exp[-jk  v  AY     ]exp[-jk    (1-u  -v  )   /ZAZ] 

=_(M_-,)/2     n    n    TM     n     o      '^       ^     J  n  o     qn       ^     J  n         o     o 


(M-D/2 

J  c    exp  [-jk  u  AX.    ]dv    du    exp[j27Tkf  t]  } ;  (2.29) 

,,,lwo     in  Jno     un       o      o  o 

m=-(  M-D/2 

u    +  v       <     1 
o         o    — 
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III.   IMPLEMENTATION  OF  THE  MODEL 

The  main  task  of  the  computer  simulation  model  is  to 
compute  and  file  the  electrical  output  signal  at  all  elements 
of  the  receive  planar  array  according  to  equation  (2.29)  .   The 
model  was  implemented  on  an  IBM  system/370  mainframe  computer 
using  FORTRAN.   The  program  is  computation  intensive  because 
the  evaluation  of  equation  (2.29)  involves  a  double  integral 
with  nested  summations.   Therefore,  much  attention  has  been 
given  to  program  efficiency.   This  section  discusses  the 
equations  which  were  actually  implemented  and  the  overall  logic 
of  the  program. 

A.   NUMERICAL  METHODS 

The  model  computes  the  real  output  electrical  signal  at 

all  receiver  elements  as  given  by  equation  (2.29) .   Note  that 

the  double  integral  with  respect  to  (w.r.t.)  direction 

cosines  u   and  v   behaves  like  a  frequency  response  or  trans- 
o       o  ^     J  c 

fer  function  H(f,i,q) ,  which  depends  on  the  position  of  the 
receiver  elements  as  indicated  by  the  indices  (i.q) .   We  define 


A   -2  1    1    (N-D/2 

H(f,i,q)   =  -2  /    /       I          *VVVn^  • 

c  -1   a    n=- (N-D/2 

o  q 


2     2  1/2  (M-1)/2 

exp[-jk  v  AY     ]exp[-jk    (1-u  -v  )    '    AZ]  )  c  exp[-jk  AX.   u  ]dv  du 

^     J  n  o     qn       ^     J  n         oo  ,£„  , ,  /0     m     ^     J  n     im  o       oc 

^  m=-(M-l)/2 

(3.1) 
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2 

Note  that  the  factor  1/c   in  equation  (2.29)  was  moved  in 

2 
front  and  replaced  by  the  factor  1/c  .   This  is  possible 

because  the  differences  in  sound  speed  at  the  different  trans- 
mitter elements  are  very  small  for  all  practical  cases  and 
the  effect  on  the  amplitude  of  H(f,i,q)  can  be  neglected. 
However,  the  dependence  of  the  initial  speed  of  sound  on  the 
index  n  of  the  transmit  array  will  be  taken  into  account  in 
the  evaluation  of  the  phase  terms  as  expressed  by  the  different 
values  of  k   in  equation  (3.1)  . 

Using  the  definition  of  H(f,i,q),  the  expression  for  the 
real  output  electrical  signal  at  a  receiver  element  with 
indices  (i,q)  becomes 


KMAX 
y(t,i,q)   =   2  Re  {   £    cvH  (kf  ,  i  ,q)  exp  (  j  2iTkf^t)  }  (3.2) 

k  =  l 


The  computer  program  model  first  calculates  H(f,i,q) 
as  given  in  equation  (3.1) .   Then,  with  the  use  of  H(f ,i,q) , 
the  output  electrical  signal  is  readily  computed  according 
to  equation  (3.2). 

The  ocean  medium  transfer  function  H   is  possibly  random 
and  therefore  difficult  to  integrate  numerically  so  a  change 
in  the  order  of  integration  in  equation  (3.1)  is  made  to 
achieve  a  more  suitable  form  for  implementation: 


25 


2  1        (N-U/2 

H(f,i,q)      =-5-       /  J  d  H    (k     V     n,q)exp[-jk  voAY      ]    • 

c     a         n=-(N-l)/2  ^ 

o     q 

1        (M-D/2  ,/? 

f  V  c  exp[-jk  u  AX.    ]exp[-jk    (1-u  -v  )    7    AZ]du     dv       (3.3) 

J  ./,,  /n     m    v     J  n  o     un      *     J  n         o     o  o       o 

-1       m=-(M-l)/2 


In  order   to    improve   efficiency,    we    concentrate    first   on 

the   inner    integral   w.r.t.    u    ,    which   we    rewrite    as 
^  o 


/      Sv(f,u    )exp{-k    [u    AX.    +    ( l-u2-v2) 1/2AZ] }du  (3.4) 

_1J         X'o^noi  oo  o 


where 


AX.       =      x      -    x      +    id' 
i  r  o  x 


Note   that 


(M-l)/2 
Sv(f/u    )       =  y  c    exp[jk   u   md    ]  (3.5 

X  °  m=-(M-D/2      ra  n   °      x 


is    the   transmit    far-field   directivity    function   of    a    line   array 
consisting   of   equally    spaced   point-sources    as    a    function   of 
direction   cosine    u      and    frequency    f.       Now   equation    (3.3)     reduces 
to: 


2  1        (N-D/2 

H(f,i,q)      =     -j       /  I         d^k     v     nfq,exp[-jknvoAY     ]    • 

c0  a         n=- (N-D/2  M 

/     S    (f,u  )exp{-jk    [u  AX.    +   (l-u2-v2) 1/2AZ] }du     dv  (3.6) 

_.;        x'       o       ^     J  n     o     l  o     o  o       o 
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In  the  case  of  simple  amplitude  weighting  in  the  x-direction 
such  as  rectangular,  triangular  or  Hamming  weighting,  we  can 
obtain  closed  form  expressions  for  S  (f ,u  ) .   This  improves 

A      O 

program  efficiency.   For  example,  in  the  case  of  rectangular 
amplitude  weighting 


sin[k  (u  -u,  )Md  /2] 
<r     t  j=         \  n   o    b    x  /-,-, 

S„(f,u  )   =   — i — n — ? n — r^n —  (3.7 

X   '  o       sin[k  (u  -  u,  d  /2] 

n   o    b   x 


where  u,  is  the  direction  cosine  value  to  which  the  transmit 
b 

far-field  beam  pattern  is  steered. 

Before  presenting  more  details  on  how  H(f,i,q)  is  calcu- 
lated, a  brief  overview  on  numerical  and  analytical  integration 
techniques  which  were  used  is  given. 
1 .   Integration  Techniques 

a.   Method  of  Stationary  Phase 

This  method  approximates  an  integral  of  the  form 


b 
I   =    /   S(z)  exp[jg(z)]dz  (3.8) 

a 


where  g(z)  is  a  rapidly  varying  phase  term  compared  to  the 
slowly  varying  amplitude  term  S (z) .   Officer  [Ref.  10]  shows 
that  this  integral  can  be  approximated  by  using  the  fact  that 

the  main  contribution  exists  at  the  stationary  point  z    which 

1    *  sp 

is  defined  by 


dg(Z   ) 

P    =0  (39) 

dz        u  ^'^' 
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The  limits  of  the  integral  may  be  extended  to  infinity 
because  only  the  region  around  the  stationary  point  contributes 
to  the  integral.   Assume  that  around  the  stationary  point,  z 
is  given  by 


z   =   z    +  c  (3.10) 

sp 


The  approximation  of  the  integral  is  then  given  by 


I   =   S(zsp)  /27T/|g"(zsp)  |exp{j[g(zsp)±^]}         (3.11) 


under  the  condition 


sg'"  (z  )  | 

^?r—   <<   1  (3.12) 


3g" (z  ) 

sp 


where  the  sign  in  the  phase  term  of  equation  (3.11)  equals 

the  sign  of  the  second  derivative  g" (z   ). 

^    sp 

b.   Numerical  Integration  Using  Newton-Cote's  Formulas 
This  standard  technique  [Ref.  11]  is  the  method 
used  by  the  program  to  perform  a  numerical  integration  over 
a  subinterval  of  the  total  region  of  integration.   Subinter- 
vals  can  be  defined  by  dividing  the  total  region  of  integra- 
tion into  equal  parts  or  by  a  scheme  as  described  in  the  next 
section.   Each  subinterval  is  integrated  by  evaluating  the 
Simpson's  rule  estimate: 


S   =   -=-(s   +  2s,  +  4s„  +  2s_,  +  4s.  +  ...  +  2s   ,  +  s 

Jo      1      2      3      4  n-1     n 


(3.13) 
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where  s  ,...,s   denote  values  of  the  integrand  evaluated  at 
on  r 

different  points  of  the  subinterval,  and  the  distance  between 
these  points  is  given  by  the  step  size  h.   The  original  step 
size  h  is  equal  to  half  the  size  of  the  subinterval  to  be 
integrated.   The  step  size  is  then  halved  to  obtain  a  more 
accurate  estimate  which  takes  more  points  into  account.   This 
is  repeated  until  a  tolerance  criterion  is  met.   Using  succes- 
sive Simpson's  rule  estimates,  one  can  do  a  Richardson's 
extrapolation  [Ref.  12].   This  extrapolation  uses  the  last 
two  estimates  by  Simpson's  rule,  S~  and  S-,  ,  where  S~  uses 

half  the  step  size  of  S, .   Both  estimates  have  a  global  error 

4 
of  order  h  .   The  extrapolated  value  becomes : 

Extrapolated  value   =   so  +  (S2-S,)/15       (3.14) 

with  an  error  of  order  h  .   If  two  successive  extrapolated 
values  are  within  tolerance,  the  last  extrapolated  value  is 
taken  as  the  value  of  the  integral  for  that  subinterval.   This 
scheme  creates  a  semi-adaptive  integration  procedure  because 
different  subintervals  can  require  different  step  sizes.   It 
could  be  further  enhanced  but  would  then  require  more  overhead 
Another  possibility  would  be  the  use  of  a  Gaussian  quadrature 
rule  which  needs  less  function  evaluations  for  the  same  order 
of  error  but  has  the  disadvantage  that  the  points  are  not 
identically  spaced  and  the  scheme  cannot  be  made  adaptive  in 
a  computation  efficient  manner. 
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c.   Numerical  Integration  Using  the  Euler  Summation 
A  problem  in  evaluating  an  integral  like  equation 

(3.4)  is  convergence.   The  real  and  imaginary  parts  of  the 
integrand  of  equation  (3.4)  are  oscillatory.   Successive 
positive  and  negative  contributions  tend  to  cancel,  but  each 
contribution  itself  is  not  negligible  and  the  integral  behaves 
oscillatory  as  a  function  of  its  limits  of  integration.   There- 
fore, one  has  to  integrate  over  a  relatively  large  region  to 
obtain  convergence.   This  is  in  conflict  with  the  desired 
computational  efficiency. 

Squire   [Ref.  13:  pp.  172-173]  describes  an  inte- 
gration procedure  for  oscillatory  integrands  which  splits  the 
range  of  integration  into  parts  using  the  zero  crossings  of 
the  integrand  as  points  of  separation.   The  individual  subin- 
tervals  can  be  evaluated  with  the. use  of  standard  methods 

(as  described  in  III.A.l.b)  and  summed  to  form  the  value  of 
the  integral.   A  technique  like  the  Euler  summation  [Ref.  13: 
pp.  172-173]  can  be  used  to  improve  convergence  of  this  sum- 
mation.  The  Euler  summation  obtains  from  an  original  sequence 


al'a2,a3' *  *  * ,an 


a  second  sequence 


bl   =   a2  +  al'    b2   =   b3+a2'  ••"  bk  =  ak  +ak+1 


and  a  third  sequence 
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Cl      =      b2    +    bl'  C2       =      b3+b2'     ••"    Ck      =      bk+l+bk 


and    so   on.       It    can   be    shown    that    the    sum   of 


11.1  1 

2  al     +    4bl    +    8  Cl    +    •'•    +   ^nml 


is  the  same  as  the  sum  of  the  original  sequence.   Convergence 
however  is  much  improved  which  cuts  down  on  the  number  of 
subintervals  to  be  integrated  and  improves  efficiency. 

This  method  works  well  on  integrals  where  successive 
contributions  from  subintervals  to  the  integral  are  of  decreas- 
ing magnitude  and  of  alternating  sign  as  illustrated  in  Figure 
2.   To  use  this  method  it  must  be  possible  to  solve  for  the 
zero  crossings  of  the  integrand,  which  generally  requires  a 
closed  form  expression  to  solve  for  the  roots  of  the  integrand. 
2 .   Application  of  Integration  Methods 

The  methods  introduced  in  the  previous  section  are 
used  to  define  a  procedure  to  calculate  H(f,i,q).   First, 

the  integration  w.r.t.  direction  cosine  u   is  considered, 

o 

which  is  given  by 


/   Sx(f,uo)exp{-jkn[uoAX.+(l-u2-v2)1/2AZ]}duo      (3.4 


Both  the  method  of  stationary  phase  and  the  numerical 
method  (Sections  III.A.l.b  and  III.A.l.c)  using  the  Euler 
summation  are  implemented.   To  apply  the  method  of  stationary 
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INTEGRAND 
G(  Z  ) 


2nd  SUBINTERVPL 


1st  SUBINTERVGL 


Fig.  2    Oscillatory  Integrand  Showing  Decreasing 
Contributions  to  the  Integral 
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phase   we    define    in    accordance   with    equation    (3.8) 


g(u    )       =      -k     [AX.u     +  (1  -u2  -  v2)  1/2AZ]  (3.15) 

3      o  n        i   o  o        o 


The  amplitude  term  S  (f,u  )  is  the  directivity  function  which 

e  x     o  2 

is  indeed  slowly  varying  w.r.t.  g(u  ) .   We  find  the  sta- 
tionary point  UOSP  by  setting  g'(UOSP)  equal  to  zero: 


g'  (UOSP)   =   -k  AX.  +k  AZ(1  -u2  -  v2)  1/2u    =   0   (3.16) 

n   l    n        o    o       o 


leading  to 

AX.      1  -v2     1/2 

UOSP   =  -~{ - j)  (3.17) 

116     1  +  (AXi/AZ) 

For  the  values  of  the  second  and  third  derivatives  one  finds 


g"(u  )   =   k  AZ(l-u2-v2)  3/2(l-v2)         (3.1! 
J    o       n       o    o  o 


g"'  (u)  =   3k  AZu  (1  -u2  -v2)  5/2(l-v2)     (3.19) 
J     o       n    o     o    o  o 


We  now  form  the  ratio  of  the  third  and  second  derivatives 
at  u   equal  to  the  stationary  point  UOSP: 

2     2  1/2 
a"'  (UOSP)               2  -l/2^AXi  +AZ  ' 
V  (will       "   3AX.(l-v2,  I/2         2 (3.20) 

A  ^j 


The  condition  in  equation  (3.12)  becomes 
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1  +  (AX./AZ) 


<<  i 


(3.21) 


from  which  we  conclude  that  the  above  condition  is  obeyed 

only  when  AZ  >>  AX.,  which  means  that  the  method  is  only 

valid  for  relatively  small  cross  ranges  x  -x   since 

■*  ^     r    o 

AX.  =  x  -x  +  i  d1  . 

i     r    o     x 

The  program  implements  the  following  procedure  to 
compute  the  magnitude  and  phase  of  the  integral  w.r.t.  u 
[see  equation  (3.4)],  denoted  by  IWRTUO. 
1)   Calculate  the  stationary  point  UOSP 


UOSP   = 


AX±   /     1  -v 
AZ  V  1  +  (AX./AZ)2 


3.17) 


2)   Calculate  the  value  of  the  second  derivative  DRV2ND 
at  the  stationary  point: 


DRV2ND   = 


2    2  -3/2      2 
k  AZ(1  -UOSP  -  v  )  J/*(l  -v  ) 
o  o  o 


(3.22) 


3)   Calculate  magnitude  and  phase  of  the  integral 


MAG   = 


S  (f,UOSP)  /2^/DRV2ND 


(3.23) 


PHASE   = 


-k  [AX. UOSP  +  (1  -UOSP2  -v2)1/2AZ]  +£ 
n    l  o  4 


Note  that  the  integral  is  a  function  of  v   and  the  indices  i 

r  o 

and  n  through  k   and  AX..   The  dependence  of  the  magnitude 
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on  the  index  n  is  negligible,  thus  DRV2ND  is  calculated  for 

all  n  using  k  .   We  will  denote  the  integral  by  IWRTUO(v  ,k  ,i) 
3   o  3      J  on 

The  numerical  method  using  the  Euler  summation  can 
also  be  applied  to  the  integration  w.r.t.  u  .   This  numerical 
method  is  implemented  to  verify  the  result  of  the  stationary 
phase  method  and  extend  application  of  the  program  to  larger 
cross  ranges.   The  applicability  of  the  method  can  be  shown 
by  separating  the  integrand  of  equation  (3.4)  into  a  real 
and  imaginary  part 


/   (S  (f ,u  )cos{-k  [u  AX.  +  (1  -u2  -  v2)  1/2AZ]  } 
,;     x    o       n   o   i       o    o 

(3.25) 

+  jS  (f  ,u  )sin{-k  [u  AX.  +  (1  -u2  - v2) 1/2AZ]  } ) du 
Jx'o       noi       o    o  o 


As  shown  above,  the  point  of  main  contribution  to  the 
integral  in  equation  (3.25)  is  the  stationary  point  UOSP. 
Away  from  this  point  the  integrand  becomes  more  and  more 
oscillatory  and  consequently  more  difficult  to  integrate 
numerically.   One  would  therefore  like  to  use  a  scheme  in 
which  the  size  of  the  subintervals  is  dependent  on  the  behavior 
of  the  integrand.   The  method  described  in  III.B.l.c  uses 
the  zero  crossings  of  the  integrand  as  points  of  separation, 
thus  making  sure  that  the  integrand  within  a  subinterval  is 
"well-behaved"  and  can  be  efficiently  integrated  by  standard 
numerical  methods.   This  division  into  "well-behaved"  subin- 
tervals could  also  be  obtained  by  using  the  relative  extremes 
as  points  of  separation.   The  integrand  of  equation  (3.25) 
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has  closed  form  solutions  to  the  zero  crossings  and  relative 
extremes  of  both  the  real  and  imaginary  parts.   The  zero 
crossings  of  the  imaginary  part  of  the  integrand  are  found 
by  solving: 


2    2  1/2 
-k  AX.u   -  k  AZ(1  -u  -v  )  '         =      nTT  ;   n  integer    (3.26) 
n   i  o     n        o    o  3 


or 

-AX.(mr)  ±  AZ  [  (1  -  v2)k2(AX2  +  AZ2)  -(n^)2]1/2 

l               o   n    l  .  _  __ . 

u    =   ~ ~ (3.27) 

°  k  (AXZ  + AZ  ) 

n    i 

These  also  are  the  locations  of  the  relative  extremes 
of  the  cosine  term  appearing  in  the  real  part  of  the  integrand. 
In  order  to  avoid  integrating  the  real  and  imaginary  parts 
of  the  integral  separately,  the  same  subintervals  are  used 
for  both  the  real  and  imaginary  parts  of  the  integral.   Points 
of  separation  for  the  subintervals  are  the  zero  crossings  of 
the  imaginary  part  of  the  integrand,  which  coincide  with  the 
relative  extremes  for  the  real  part  of  the  integrand.   As  a 
starting  point  for  the  integration,  define  the  first  subinter- 
val  using  the  point  of  main  contribution,  i.e.,  the  stationary 
point  UOSP. 

Fig.  3  shows  a  typical  behavior  of  the  real  and  imagin- 
ary parts  of  the  integrand  around  the  stationary  point.   The 
figure  also  indicates  the  position  of  the  subintervals. 
Integration  of  the  successive  subintervals  (from  UOSP  outward), 
will  form  alternating  terms  of  decreasing  magnitude  for  both 
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IMAGINARY 


REAL 


PHASE 


Fig.    3         Typical   Behavior   of    Integrand    in 

Eq.     (3.25)    Around   Stationary   Point 
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the  real  and  imaginary  parts  of  the  integral,  making  it 
possible  to  effectively  speed  up  the  convergence  by  using  a 
complex  Euler  summation.   The  decreasing  magnitude  is  caused 
by  the  non-linear  behavior  of  the  phase  term  and  the  decreas- 
ing amplitude  of  the  integrand  which  is  given  by  the  transmit 
far-field  directivity  function  S  (f  ,u  ) .   The  program  imple- 
ments the  following  procedure: 

1)  Calculate  the  value  of  the  stationary  point  UOSP 
according  to  equation  (3.17).   This  value  of  u   also 
denotes  the  maximum  (negative)  value  of  the  phase 
term. 

2)  Calculate  NZ  the  maximum  integer  multiples  of  tt 
contained  in  the  maximum  (negative)  value  of  the 
phase  term  (see  Fig.  3) . 


NZ       =       INTEGER ( 


-k     [UOSPAX.  +  (l-UOSP2-v2)  1//2AZJ 
n i o 

TT 


)   (3.28) 


Note  that  NZ  denotes  the  zero  crossings  of  the  imaginary 
part  of  the  integrand  nearest  to  the  stationary  point 
UOSP. 

3)   Calculate  the  values  of  u0  for  the  zero  crossings  of 
the  imaginary  part  of  the  integrand  denoted  by  NZ 
using  equation  (3.27)  . 


UOZ 


•AX.  (NZtt)  ±AZ[  (l-v2)k2(AX2+AZ2)-(NZTr)  2]1/2 
l  o   n    l 


1/2 


k  (AX2  +AZ2) 
n    i 


(3.29) 


Integrate  the  subinterval  between  these  two  values 
of  UOZ. 


4)   Decrement  NZ  and  calculate  the  values  of  UOZ 


1,2 


denoting  the  next  zero  crossings  of  the  imaginary 
part  of  the  integrand,  by  using  the  formula  given 
in  3)  above. 


38 


5)  Define  two  subintervals ,  one  on  each  side  of  the 
stationary  point,  from  zero  crossing  to  zero  crossing 
for  the  imaginary  part  and  from  extreme  to  extreme 
for  the  real  part  of  the  integrand. 

6)  Integrate  both  subintervals  and  add  the  contribution 
to  the  total  value  of  the  integral  using  a  complex 
Euler  summation. 

7)  Check  for  convergence  by  comparing  the  difference 
between  the  last  two  values  of  the  Euler  summation 
with  a  tolerance.   If  not  convergent,  go  back  to  4) 
and  repeat  the  procedure. 

One  more  note  can  be  made  concerning  the  numerically 

calculated  value  of  the  integral  w.r.t.  u   for  different  k  . 

J  o  n 

Recalculation  of  the  whole  integral  for  different  values  of 
n  would  be  very  time  consuming.   We  therefore  introduce  a 
correction  term  to  obtain  values  of  the  integral  for  n  ^  0 , 
based  on  the  value  of  integral  for  n  =  0.   This  correction 
term  is  only  of  importance  to  the  phase  of  the  integral,  the 

effect  of  different  values  for  k   on  the  magnitude  is 

n  3 

neglected.   For  n  =  0  the  stationary  phase  method  results  in 
the  following  expression  for  the  phase: 


PHASE  (n  =  0)   =   -ko[AX.UOSP+(l-UOSP2-Vo)  1//2AZ]+  j       (3.30) 


Brekhovskikh  [Ref.  14]  demonstrates  that  taking  into 
account  higher-order  terms  in  the  method  of  stationary  phase 
does  change  the  amplitude  but  not  the  phase  of  the  resulting 
approximation.  This  shows  that  the  point  of  main  contribu- 
tion, i.e.,  the  stationary  point  determines  the  phase  of  the 
integral.   If  the  method  of  stationary  phase  is  no  longer 
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valid  we    can  write    a    general    form    for    the   phase   based    on 
equation    (  3 . 30 )  : 


PHASE     (n=0)       =      -k     [AX.UOSP+(l-UOSP2-v2)  1//2AZ]     +    £   +    A<j> 

O    1  O  4 


(3.31) 


where  Ac{)  represents  some  unknown  correction  term.   This  relies 
on  the  fact  that  the  stationary  point  UOSP  is  always  very 
close  to  the  point  of  main  contribution.   Based  on  this,  the 
phase  term  for  n  jt   0  can  be  written  as: 


PHASE(n)  =  -[k  +(k  -k  )  ]  [  AX  .  UOSP+  ( l-UOSP2-v2)  1/2AZ]  +x  +  A<j>   (3.32) 
o    n   o      l  o        4 


where 


(k  -k  ) [AX. UOSP  +  (l-UOSP2-v2)1/2AZ]  (3.33) 

n   o     i  o 


represents  the  correction  term  which  is  readily  computed  for 
different  values  of  n.   This  phase  correction  term  has  to 
be  added  to  the  phase  for  n  =  0  to  obtain  phase  values  for 
n  ^  0.   Note  that  this  correction  term  is  only  important  for 
larger  ranges  AZ,  assuming  that  AZ  >>  AX.   This  can  be  seen 
by  calculating  a  typical  value  for  k  -k  ,  where  the  element 
spacing  is  taken  as  A/2  and  the  sound  speed  profile  has  a 
constant  gradient  g: 


k  -k    =   2TTf(  — — )  z      -  ^2  n   =   3.6  x  10  5  n         (3.34) 

no  c     c  c 

no  o 
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where  we  used  c   =  1500  m/sec  and  g  =  0.017  1/sec.   An  upper 

limit  for  the  correction  term  is  (k  -k  )AZ  if  AZ  >>  AX. 

n   o 

Having  established  two  methods  to  implement  the  calcu- 
lation of  the  inner  integral  w.r.t.  direction  cosine  u  ,  the 

outer  integral  w.r.t.  direction  cosine  v   is  now  considered. 
^  o 

The  integral  can  be  written  as  [see  equation  (3.  6) J  : 


1    (N-D/2 
/        y      d  H  (k  ,v  ,n,q)exp[-jk  v  AY   ]IWRTUO(v  ,k  , i ) dv 
a    n— (N-D/2   n  M  n   °  n  °   qn         °   n     ° 

q 

(3.35) 

The  ocean  medium  transfer  function  values  are  calculated  with 
the  use  of  equations  (2.11) ,  (2.16)  and  (2.18) .   Note  that  the 
deterministic  and  random  phase  terms  pertain  to  the  case  of 
a  sound  speed  profile  with  a  constant  gradient.   In  applying 
these  formulae,  we  have  to  take  into  account  that  the  source 
depth  for  a  transmit  element  is  not  y   but  y  +ndy .  This  results 
in  the  following  formulae  which  were  actually  implemented: 


AM   "   (koVo»"1/2  (3-36) 

k   c      c 
9MD(n^>   "   2^'f  (c(yr  +  qa;>  -  1}     +  AYqn'  (3-3?> 

eMR(n^»   "   !r'Tg"NR  ln<c(y  !qd')"  <3-38) 

o   ^  2  r  ^  y 


We  note  that  both  the  amplitude  and  phase  terms 

depend  on  the  value  of  the  wavenumber  k   and  direction  cosine 
r  n 
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v  .   The  expressions  within  the  square  brackets  for  9..^  and 
o  MD 

9    are  only  dependent  on  the  geometry  and  can  be  precomputed. 
A  random  number  generator  was  used  to  generate  a  series  of 
N(0,1)  distributed  numbers  to  be  used  for  n    in  equation 
(3.37) . 

Because  of  the  complexity  of  the  total  phase  term  in 
the  integrand  of  equation  (3.35),  a  closed  form  approximation 
could  not  be  found.   The  numerical  method  used  is  also  based 
on  the  procedure  given  in  III. A. I.e.   The  integral  is  divided 
into  subintegrals  in  such  a  way  that  the  subintegrals  are 
all  "well  behaved."   The  problem  is  to  find  a  general  form 
which  approaches  the  phase  of  the  integrand.   This  form  can 
then  be  used  to  calculate  the  points  of  separation  in  order 
to  define  the  subintervals .   We  will  proceed  by  separating 
the  phase  term  of  the  integrand  into  a  major  part  and  a  minor 
slowly  varying  part.   Combining  the  phase  term  of  the  integral 
w.r.t.  u   given  by  equation  (3.32)  with  equation  (3.35) 
results  in  the  following  expression  for  the  integration  w.r.t. 

vo: 


1       (N-D/2 
/  d  H.,(k   ,v   ,n,q)exp{-jk  v    (qd*  -nd  )  +  T  +  Acj)} 

a;     n=-(N-l)/2     n^-    n     °  n  °       *         *       4 

exp{-j(k  -k   )  [AYv  +AX.UOSP+(l-UOSP2-v2) 1^2AZ] }    • 
^     J     n     o  o       l  o 


IWRTUO(v   ,k   ,i)  |exp{-jk    [v  AY+UCSPAX .  +  ( 1-UOSP  -v2)AZ]}dv  (3.39 

o    n       '     *    J  o    o  i  o  o 
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where    AY    =    y     -  y    .       If   we    let 
2  r      2  o 


(N-D/2 

X(v  )      =  V  d  H__(k    ,v   ,n,q)exp{-jk  v   (qd'  -nd  )  +  7-  +    Ad)}- 

0  n=-(N-l)/2     n  ™     n     °  Jno^y         y'      4 


2      2   1/2 

exp{-j(k  -k   )  [AYv  +AX . UOSP+ ( 1-UOSP  -v  )    x    AZ]}-    IWPTUO(v   ,k    ,i)  (3.40 

c    J     n    o  o       1  o  '  on' 


and 


a(v    )       =      -k     [AYv     +UOSPAX.    +  ( l-UOSP2-v2) 1/2AZ] 
0001  o 

(3.41 

a(v    )       =      -k    [AYv     +(l-v2)1/2   Ax2  +AZ2] 
o  000  1 


Then    the    integral    given   by   equation    (3.39)    becomes 


1 
/      X(vq) exp[ ja(vQ) ]    dvQ 
a 

q 

where  X(v  )  is  slowly  varying  compared  to  exp{ja(v  )}.   The 

phase  term  a(v  )  can  now  be  used  to  define  subintegrals  which 
o  J 

are  "well  behaved."   This  is  accomplished  by  taking  the  roots 

of  sin(a(v  ))  or  cos(a(v  ))  as  the  points  of  separation  between 
o  o  c  i. 

the  different  subintervals .   Contributions  from  the  subintegrals 

are  in  general  of  decreasing  magnitude  and  alternating  sign, 

so  the  Euler  summation  can  be  used  to  effectively  speed  up 

convergence . 

The  implemented  procedure  uses  the  roots  of  sin(a(v  )) 
c  c  o 

as  the  points  of  separation  and  is  similar  to  the  one  already 
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described  for  the  integration  w.r.t.  u  .   The  point  of  main 

3  o        c 

contribution  is  approximately  where  the  phase  term  a (v  )  is 
stationary.   This  point  is  used  to  begin  the  numerical  inte- 
gration.  The  stationary  point  VOSP  is  given  by 


VOSP   =   AY/(R2.  +  AY2)1/2  (3.43) 


which  is  a  solution  of  a1 (VOSP)  =  0  where  R  .  is  the  radial 

ri 

distance  from  the  center  of  the  transmit  array  to  a  receiver 
element  with  index  i  and  is  given  by 


2      2  ^2 
R  .   =   (AX   +  AZ  )  (3.44) 

ri        l 


The  points  of  separation  for  the  integral  can  be  found  by 

1/2 
-AY(n7r)±R  .  [k2(AY2  +  R2.  )  -(nir)2] 
VOZ      =   E± 2 £i_^ (3.45) 

ko(Rri  +  AY  } 

Another  method  using  a  brute  force  technique  was  used 
to  validate  the  above  described  procedure.   This  method 
integrates  w.r.t.  direction  cosine  v   between  user  specified 
limits  of  integration.   It  divides  the  range  of  integration 
into  subintervals  of  equal  size  which  are  integrated  by  a 
standard  routine  as  described  above  in  III.A.l.b.   Both 
methods  agreed  numerically  although  the  method  using  the 
unequally  spaced  subintervals  and  the  Euler  summation  was 
much  more  efficient  in  terms  of  computation  time. 
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B.   PROGRAM  OUTLINE 

The  flowchart  in  Fig.  4  gives  an  overview  of  the  program 
logic.   This  logic  is  controlled  by  flags  which  are  read 
from  the  input  data  file  together  with  all  program  parameters 
The  program  parameters  determine: 

-   Model  geometry:   x  ,  y  ,  z  ,  x  ,  y  ,  z  ,  M,  M,  M '  ,  N1, 

d  ,  d  ,  d',  d' .    °    °    °    r    r    r 

x    y    x    y 

Speed  of  sound  profiles:  c    ,    g ,    a . 

Beam  steering  angles  for  the  transmit  array  and  the 
frequency  to  be  used  for  the  beam  steering  calculations. 

The  fundamental  period  TQ  of  the  electrical  input  signal 
and  the  total  number  of  samples  L  in  the  electrical 
output  signal.   The  program  determines  the  sampling 
period  Tg  =  TQ/L  and  calculates  the  output  electrical 
signal  at  time  instants 


-(L-l)T  /2,  ...,  0,  ...,  (L-l)T  /2 
o  o 


The  values  of  the  frequency  components  in  the  electrical 
input  signal  and  their  complex  Fourier  series  coefficients 

V 

Tolerance  values  ERRORU  and  ERRORV  for  the  numerical 
integrations  w.r.t.  u   and  v  .   The  user  should  validate 
his  choice  for  the  tolerance  values  by  comparing  them 
against  the  absolute  values  calculated  by  the  program 
for  the  integrations. 

Processing  is  done  in  accordance  with  flag  settings  and 

allows  for  a  number  of  options: 

-   The  frequency  response  H(f,i,q)  can  be  calculated  or 
read  from  a  previously  created  file. 

The  calculated  frequency  response  H(f,i,q)  can  be 
filed  for  later  use  and/or  printed. 

The  output  electrical  signal  can  be  calculated  and 
filed  for  further  processing  by  other  programs. 
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(    st'rt    ) 


Calculate 
H(f,  i,q) 


Read  input 
from  file 


(^    Stop ) 


f     Stop     J 


/Read 
H(f,  i,q> 
from  file 


Print 
H(f,  i,q) 


'Write 
H(fi,  ,q) 
to  file 


Calculate 
output  signal 
Y(t) 


Write 
Y(-£,  i,q) 

to  f i le  i 


Fig.  4    Flowchart  of  Main  Program  Module 
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Most  important  numerical  calculations  required  by  the 

model  are  done  within  the  module  CALCHF.   The  module  CALCHF 

calculates  the  frequency  response  matrix  H(f,i,q),  i.e.,  the 

frequency  response  for  all  receiver  elements  (i,q).   This 

mainly  involves  evaluation  of  the  double  integral  w.r.t.  u 

o 

and  v  .  It  has  two  versions  dependent  on  the  integration 
technique  used  to  evaluate  the  integral  w.r.t.  u  : 

CALCHF__1    Use  stationary  phase  method  to  integrate 

w.r.t.  u  . 
o 

CALCHF_2    Use  numerical  technique  to  integrate  w.r.t.  u  . 

To  avoid  recalculation  and  improve  program  efficiency 

the  outer  integration  w.r.t.  v   is  done  simultaneously  for 

o  J 

all  elements  with  the  same  index  i.   This  is  particularly 
important  when  using  the  more  time  consuming  numerical  method 
to  integrate  w.r.t.  u  .   This  can  be  done  by  separating  the 
summation  inside  the  integrand  into  a  part  dependent  on  q 
and  a  part  dependent  on  i.   Recall  that  equation  (3.35) 
gives  the  form  of  the  integrand  as  follows: 


(N-l)  /2 

V     dR  (k  ,v  ,n,q)exp{-jk  v  AY  }lWRTUO(v  ,k  ,i)         (3.46) 
n=-(N-l)/2  n^  n  °  n°  ^        °  n 


So  given  a  value  for  IWRTUO  for  a  given  index  i  we  can  evalu- 
ate the  integrand  for  all  values  of  q,  thus  avoiding  recalcu- 
lating the  most  computation  intensive  part.  This  scheme  does 
result  in  more  complex  coding.   The  program  has  to  do  N1 
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integrations  and  keep  track  of  N'  Euler  summations  simul- 
taneously.  It  does,  however,  also  reduce  some  overhead. 
Convergence  checking  on  the  integration  w.r.t.  v   is  done 
only  on  the  two  outermost  elements  with  indices  q  =  (N'-l)/2 
and  q  =  -(N'-l)/2. 

The  calculation  of  those  terms  that  multiply  IWRTUO(k  ,v  ,i) 

no 

in  the  summation  given  by  equation  (3.46)  is  simplified  by 
precomputing  part  of  the  phase  components  of  the  ocean  medium 
transfer  function  given  by  equations  (3.36)  and  (3.37).   The 
program  does  this  in  the  subroutine  HMWKB  and  stores  the 
resulting  phase  integral  values  in  an  array.   To  obtain  a 
particular  value  for  the  phase,  the  program  only  needs  to 
multiply  with  the  factor  k  /v  .   Especially  in  the  case  of  a 
more  complicated  sound  speed  profile,  the  computation  of 
the  phase  integrals  can  be  involved  and  precomputation  im- 
proves efficiency.   The  stated  time-invariance  of  the  model 
also  forces  precomputation  since  the  random  numbers  used  to 

model  the  phase  term  0.,_  can  be  drawn  only  once.   The  random 

c  MR  ■* 

numbers  were  obtained  from  the  IMSL  routine  GNNML  which  was 
used  to  generate  a  series  of  N  xN1  random  numbers  with  distri- 
bution N(0,1).   From  these  random  numbers  the  different  random 
phase  terms  9    for  all  combinations  (n,q)  are  readily  com- 
puted  using  equation  (3.37) .   Figure  5  shows  a  simple  flow- 
chart of  the  CALCHF  module. 

The  innermost  block  which  actually  calculates  H(f,i,q) 
is  shown  in  more  detail  in  Fig.  6.   The  diagram  shows  how  the 
subintervals  are  defined  by  the  program.   The  procedure  followed 
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CftLCHF 


Call  HMWKB 
Calculate 
phase  integral 


DO 

'for  all  frequency^ 
comDonents 


DO 

for  all  i 


Calcu 

late 

H(f, 

i,q) 

for  a 

11  a 

c 


Return 


Fiq.  5    Flowchart  of  Module  CALCHF 
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? 


Calculate 
V0SP,  NZ 


Set: 
V0BEG1  =  V0SP 
V0BEG2  =  V0SP 


Using  NZ  calculate 
zero  crossings 
V0END1, V0END2 
(Eq.  3.45) 


Define  intervals: 
right:  V0BEG1  -  V0END1 
left:   V0END2  -  V0BEG2 


Integrate 
subintervals 


Add  contribution 
to  Euler  summation 


NZ  =  NZ  +  1 
V0BEG1  =  V0END1 
V0GEG2  =  V0END2 


V0SP:  Approximate  point  of  main  contribution  of  the  integral. 
NZ   :  Maximum  negative  integer  denoting  the  multiples  of 

contained  in  the  maximum  negative  value  of  the  phase. 
Note:  The  "Integrate  Subintervals"  block  makes  a  Subroutine 

call  to  evaluate  the  integrand  given  by  Eq.  (3.46). 


Fig 


Calculation  of  Subintervals  by  CALCHF 
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is  described  in  III. A. 2  for  the  integration  w.r.t.  v   and  is 

^  o 

similar  to  the  procedure  used  for  the  numerical  integration 

w.r.t.  u   in  Section  III. A. 2. 
o 

Implementation  of  the  integration  over  a  subinterval  using 
the  Newton  Cote's  formulas  is  straightforward  and  will  not 
be  examined  in  detail.   The  framework  for  this  part  of  the 
program  can  be  found  in  Squire  [Ref.  13:  p.  281]. 

The  Euler  summation  is  used  to  speed  up  convergence  of 
the  numerical  integration.   The  objective  is  to  form  the  sum 


fl    ^1    fl    ^1 

2     4     8  +  16 


from  the  original  sequence  a,  ,  a0  ,  a_. ,  a,  ,  ...,  a   which  are  the 
r        -i        1234         n 

contributions  from  the  subintervals  1  through  n,  where 


bk   =   ak  +  ak+l'    ck   =   bv  +  bk+l' 


The  algorithm  used  is  described  below  and  uses  an  array 

COEF  to  store  the  coefficients  a,b   -,/C   0^d   ,,  ...  etc 

n    n-1    n-2    n-3 

The  size  of  this  array  determines  how  many  contributions  can 
be  handled. 

DIMENSION  C0EF(2,SIZE) 
REAL  EULSUM, FACTOR 
INTEGER  NEW,  OLD 
C      Add  contribution  N 

COEF  (NEW,1)  =  N  th  CONTRIBUTION 
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DO  100  J  =  1,  N 
100  COEF(NEW,J+l)  =  COEF(NEW,J)  +  COEF(OLD,J) 

FACTOR  =  FACTOR  /  2 

EULSUM  =  EULSUM  +  COEF (NEW , N) /FACTOR 
C     Prepare  for  next  contribution 

SWAP(  NEW,  OLD  ) 

The  program  implements  this  algoritm  in  complex  form  to 
perform  N'  summations  simultaneously. 

More  details  on  the  actual  programming  of  the  discussed 
numerical  methods  and  flowcharts,  data  structures  used,  etc., 
is  considered  applied  programming  and  will  not  be  discussed 
in  this  thesis. 
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IV.   COMPUTER  SIMULATED  DATA 

This  section  describes  the  validation  of  the  model.   It 
presents  and  interprets  the  data  generated  by  the  computer 
model  for  a  number  of  test  cases. 

A.   PRESENTATION  OF  COMPUTER  SIMULATED  DATA 

Data  generated  by  the  computer  model  consists  of  the 
frequency  response  matrix  H(f,i,q)  and  the  output  electrical 
signal  at  each  element  y(t,i,q) .   Interpretation  of  the  data 
is  done  directly  by  examining  the  frequency  response  H(f,i,q) 
and  by  the  use  of  a  three-dimensional  DFT  program  to  calculate 
the  Fourier  transforms  of  the  computer  simulated  output  elec- 
trical signals  w.r.t.  time  and  space.   The  program  is  based 
upon  3-dimensional  DFT  beamforming  for  planar  arrays  as  des- 
cribed by  Ziomek  [Ref.  4].   First  a  brief  overview  of  the 
output  from  the  3-D  DFT  program  is  given. 

The  program  calculates  the  3-dimensional  DFT  of  the 
electrical  output  signal  from  all  receiver  elements  and  is 
given  by 


M'    N'    L' 
Yq(q,r,s)   =111    c  d  y(£,m,n)exp[-j27Tq£/L] 
m=-M'  n=-N'  £=-L' 


exp[j27Trm/(M+Mz)]exp[j2^sn/(N+NZ)  ]  (4.1) 
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where 


L:   total  number  of  time  samples  in  one  fundamental 
period  T   of  the  signal. 

L1   =   (L-l)/2 

T - :   sampling  period 

M:   total  odd  number  of  elements  in  the  x-direction 
in  the  receive  array 

M'   =   (M-D/2 

N:   total  odd  number  of  elements  in  the  y-direction 
in  the  receive  array 

N'   =   (N-D/2 

d  ,  d  :   spacing  in  x-  and  y-directions ,  respectively 
x    y 

c  ,  d  :   complex  separable  weights  for  the  array 
given  by 

c      a  exp{j0  } 
m      m   ^  J  m 

d    =   b  exp{j<J>  } 
n      n      n 

q,r,s:   binnumbers  of  the  DFT  which  are  related  to 
values  of  frequency  f,  direction  cosine  u 
and  direction  cosine  v,  respectively. 

y(£,m,n) :   output  electrical  signal  at  time  instant 

IT      at  element  (m,n)  of  the  receive  array. 

MZ ,  NZ :   Number  of  zeros  used  for  padding  to  increase 
the  resolution  of  the  direction  cosines  u 
and  v,  respectively. 

Note  that  the  symbols  used  to  characterize  this  planar  receive 
array  are  not  the  same  symbols  used  to  describe  the  receive 
array  in  the  computer  model  (II.A.l) .   Rectangular  amplitude 
weighting  without  beam  steering  is  used  for  all  the  results 
presented  in  this  chapter. 
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In  II. B. 3  we  made  the  constraint  to  represent  the  transmit 
electrical  signal  by  a  finite  Fourier  series  with  maximum 
frequency  KMAX • f  .   Thus,  to  avoid  aliasing,  we  must  sample  the 
received  signal  with  a  sample  rate  equal  to  or  greater  than 
the  Nyquist  rate: 


f    >   2  KMAX  f   ;    T    >   2  KMAX  T        (4.2 
s   —         o      o   —         s 


The  total  number  of  samples  L  to  be  taken  in  order  to  avoid 
aliasing  must  obey  the  following  inequality: 

L   >   2  KMAX  +1  (4.3) 


The  sampling  period  is  determined  by 


T    =   T  /L  (4.4 

s      o 


The  program  evaluates  the  3-dimensional  summation  in 
three  steps.   The  first  step  is  to  perform  a  DFT  w.r.t.  time 
at  each  receiver  element  (m,n) .   If  the  sound  field  incident 
upon  the  array  is  a  general  plane  wave,  then  the  received 
electrical  signal  can  be  expressed  by 


y(£,m,n)   =   g ( £T  -  t  ' 

so 


where  t   =  (u  md   +  v  nd  ) /c  denotes  the  time  delay  between 
o      o   x     o   y 

the  signals  at  the  different  elements.   The  DFT  w.r.t.  time 
becomes 
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Ys(q,m,n)   -   cmdn  J_l,  g(»Vt0)wg*         (4.5) 

where  Wj1   is  the  factor  exp{-  j  2-nqi/L}  .   Using  the  above  stated 
fact  that  the  received  signal  is  given  by  a  finite  Fourier 
series,  we  obtain 


Yr,(q,m/n)   =   c  d  Lc  exp{-j2-rrqf  (u  md  +v  nd  )/c}     (4.6) 
S  ^  m  n   q   ^ L  J   M  o   o   x    o   y  ' 


A  normalized  expression  is  defined  by 


YSN(q,m,n)   =   Yg (q ,m,n , ) / (c  dnL) 


c    exp{-j27rqf    (u   md     +v   nd    ) /c }  (4.7) 

q^J^ooxoy 


So  Y   (q,m,n)  allows  us  to  estimate  the  amplitude  of  the  fre- 

O  1 M 

quency  components  in  the  signal  at  each  element  (m,n) . 
The  next  DFT  calculation  is: 


Ys(q,r,n)   =         Yg (q ,m, n) exp{ j 27rmr/ (m+mZ ) }         (4.8) 
m=-M ' 


where  we  have  padded  with  NZ  zeros  to  increase  resolution 
Substitution  equation  (4.6)  into  equation  (4.8)  leads  to 


Yc,(q,r,n)       =      Led      exp{-j2iT  qf     v    n  d    /c }     * 
S    ^'  qn         ^      J        Moo         y 


M' 

y         c    exp{-j27rqf    u   md   /c}W.™    „  (4.9 

=-m  «       m  °    °      x  M+MZ 
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and  a  normalized  expression  can  be  written  as 


YSN(q,r,n)   =   cg  exp{  -  j  2tt  q  fQ  vq  n  d^/c}  • 

M '  M1 

c   exp{-j2Trqf  u  m  d  /c}W??„»/   J   a       (4.10 
\„,   m    r   J    -1  o  o    x     M+MZ    ^  ,  m 
m=-M  m=-M 


This  expression  is  directly  proportional  to  the  Fourier 
series  coefficients  c   and  to  the  normalized  far-field  beam 

q 

pattern  at  a  frequency  q-f   of  a  line  array  consisting  of  M 

elements  which  is  lying  along  the  X-axis.   The  planar  array 

has  N  such  line  arrays  as  indicated  by  the  index  n.   If  no 

beam  steering  is  done  (c   =  a  )  ,  then  Y„.T(q.r.n)  has  a  maximum 
3  m     m  SN  ^ 

approaching  c   at  the  bin  number  related  to  the  direction 

cosine  value  which  is  closest  to  u  .   The  DFT  bin  number  r 

o 

is  related  to  the  direction  cosine  value  by 


u   =   rA/[(M+MZ)d  ]  ;     A   =   c/qf  (4.11) 


This  value  of  u  is  an  estimate  of  u  .   Therefore,  for  each 

o 

frequency  component  of  a  plane  wave  propagating  in  the  u  ,  v 

direction,  the  expression  given  by  equation  (4.10)  can  be  used 

to  estimate  the  magnitude  of  the  frequency  component  and  the 

direction  cosine  u  .   An  estimate  of  vncan  be  obtained  in  the 

o  o 

same  way  by  evaluating 
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YSN(q,m,s)       =      c      exp{-  j  2tt   q  fQ  uq  m  dx  /c} 


N'  N1 

J  d  exp{-j  2Trqf  v  nd  /c}wfT^^7„/   7    b    (4.12) 

\, ,  n   ^   J    M  o  o    y'    N+NZ'   \T.   n 

n=-N'  7          n=-N ' 


The  above  mentioned  estimates  of  c  ,  u   and  v   are  provided 

q    o       o      ^ 

by  the  program  in  the  form  of  printouts  of  Y   (q,m,n) , 
YQKr(q,r,n)  and  Y   (q,m,s).   The  printouts  are  done  for  user 
specified  values  of  m  and  n.   No  padding  with  zeros  is  done 
for  the  transform  w.r.t.  time.   Padding  for  the  spatial  trans- 
forms is  done  to  obtain  a  user  specified  step  size  in  direction 

cosine  u  or  v.   Graphical  output  for  estimation  of  u   and  v 

r  r  o      o 

was  also  obtained.   The  complete  three-dimensional  transform 
Y~(q,r,s)  was  used  to  produce  a  three-dimensional  plot  to  aid 
in  interpretation  and  for  illustrative  purposes.   The  esti- 
mates of  u   and  v  can  also  be  used  to  calculate  estimates 
o      o 

of  the  spherical  angles  denoting  the  direction  of  the  sound 
source  as  measured  from  the  center  of  the  receive  array. 
These  angles  are  also  referred  to  as  target  bearing  and  eleva- 
tion angles.   The  estimates  are  given  by 


\p        =      tan  1(v  /u  )  (4.13) 

yo  o      o 


s 


in  1( [u2  +  v2]1/2)  (4.14) 


where  iii   and  6   are  estimates  of  the  tarqet  bearing  and  ele- 
ro       o  ^  3 

vation  angles,  respectively.   Note  that  an  error  in  the 


estimation  of  either  u0  or  vQ  or  both  affects  both  the  target 
bearing  angle  and  target  elevation  angle  estimates. 

B.   TEST  CASES 

Most  test  cases  run  by  the  model  are  kept  simple  to  aid  in 
interpretation.   All  cases  are  run  first  in  a  deterministic 
homogeneous  medium  so  that  results  could  be  easily  interpreted 
and  compared  to  known  methods  and  approximations.   All  test 
cases  are  checked  to  obey  the  criteria  required  by  the  trans- 
fer function. 

The  values  of  the  basic  parameters  such  as  source  depth 
y  ,  gradient  g  and  reference  speed  of  sound  c   were  chosen  in 
order  to  model  the  wave  propagation  between  a  source  located 
on  the  SOFAR  channel  axis  and  a  deep  receiver.   At  moderate 
latitudes   from  60  °S  to  60 °N  the  speed  of  sound  on  the  SOFAR 
axis  ranges  from  1450-14  85  m/sec  in  the  Pacific  and  from 
1450-1500  m/sec  in  the  Atlantic  [Ref.  15].   The  depth  of  the 
SOFAR  axis  is  typically  between  1000  and  1200  m  at  moderate 
latitudes  [Ref.  15].   The  gradient  below  the  SOFAR  axis  is 
given  by  g  =  0.017  1/sec  [Ref.  19] .   Values  for  the  standard 

deviation  a  of  the  random  component  of  the  index  of  refraction 

-4  -4 

range  from  0.25  x 10    [Refs.  17,18]  to  1.0  x 10    [Ref.  16]. 

These  values  were  used  in  choosing  representative  parameters 

for  the  following  test  cases  : 

1 .   Homogeneous  Cases 

As  a  baseline  test  case  we  define  the  deterministic 

homogeneous  case  HMGl : 


59 


c   =  1475  m/sec . 
o  ' 

g  =  0:   constant  speed  of  sound. 

o=0:       deterministic  medium. 

M  =  11,  N  =  11 :   a  11  by  11  element  transmit  array  with 
rectangular  amplitude  weighting  in  both  the  x 
and  y  directions. 

M'=5,N'  =5:   a  5  by  5  element  receive  array. 

x=0.0     m      x=0.0      m 
o  r 

y   =  1000  .0  m      y   =  2000  .0   m 

1  o  2  r 

z   =  0.0     m      z   =  1732.05  m 
o  r 

d   =  d   =  d1  =  d1  =  0.7375  m:   spacing  equal  to  A  .  /2 . 
x     y     x     y  r-     ^   n         mm 

8   =  30°,  ip   =  90°:   spherical  angles  denoting  the  line 

of  sight  between  the  center  of  the  transmit  array 
and  the  center  of  the  receive  array. 

u   =  0,  v   =0.5:   line  of  sight  direction  cosine  values 
r       r  ^ 

from  the  center  of  the  transmit  array  to  the  center 

of  the  receive  array. 

u,  =  0,  v,  -  0.5:   direction  cosine  values  used  to  steer 
b        b 

the  main  lobe  of  the  transmit  far-field  beam  pattern 

towards  the  center  of  the  receive  array. 

KMAX  =  1:   single  frequency  component. 

f   =  1000  Hz,  T   =  0.001  sec. 
o  o 

c,  =  a-,  =1.0:   complex  Fourier  series  coefficient  equals  1. 

L  =  3:   3  times  samples  per  fundamental  period  T  . 

Based  on  the  definition  of  test  case  HMGl  we  define 

the  following  cases: 

HMG2 :   Same  as  HMGl  except 

x   =  200  m,       z   =  1720  .464  m 
r  r 

ur  =  ub  =  0.1,    vr  =  vb  =  0 . 5 

9   =  30.66°,  ip   =  78.69°  . 

r  r 
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HMG3:   Same  as  HMG1  except 

y   =  1500  m,     z   =  866.0246  m 
J  r  r 

HMG4 :   Same  as  HMGl  except 

y   =  2500  m,     z   =  2598.074  m 
2  r  r 


HMG5:   Same  as  HMGl  except 


d  =  d   =  d*  =  d' 

x     y     x     y 


.36875  m 


f   =  2000  Hz,    T   =  0.0005  sec. 
o  o 

HMG6 :   Same  as  HMGl  except 


HMG7 


x   =  1224.745  m 
r 

u   =  u,  =0.6124, 
r     b  ' 


9   =52.24°, 
r 

Same  as  HMGl  except 


z   =  1224.745  m 
r 


v 


v,  =  0.5 
b 


\b      =    39.23° 


1-IMG8 


d   =  d   =  d'  =  d'  =  0.184375  m 
x     y     x     y 


f   =  4000  Hz, 
o  ' 


T   =  0.00025  sec. 
o 


Same  as  HMG2  except 

KMAX  =  5 :   five  frequency  components 


f   =  1000  Hz, 
o  ' 

f  =  5000  Hz. 
max 


T   =  0.001  sec. 
o 


f,  =  3000  Hz:   frequency  used  for  steering  the 

transmit  beam  pattern. 

d   =  d   =  d'  =  d'  =  1475/2  f     =  0.1475  m 
x     y     x     y  max 

c,  =  a,  =  1.0;   k  =  1,2,...,  5:   all  complex  Fourier 

series  coefficients  equal  1. 

L  =  11 :   11  time  samples  per  fundamental  period  T  . 

2 .   Deterministic  Inhomogeneous  Cases 

The  following  deterministic  inhomogeneous  cases 
incorporate  a  sound  speed  profile  with  a  constant  gradient. 
The  definitions  are  based  on  those  used  for  the  homogeneous  cases 
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INHMG1,  INHMG2,  and  INHMG5 :   Same  as  HMGl ,  HMG2  and  HMG5  except 

g  =  0  .017.  1/sec. 
INHMG8:   Same  as  HMG8  except 

g  =  0.017  1/sec. 

3 .   Random  Inhomogeneous  Cases 

The  following  cases  with  an  index  of  refraction  com- 
posed of  a  deterministic  and  a  random  component  are  defined: 
RINHMGl,  RINHMG2 :   Same  as  INHMGl ,  INHMG2  except 

a  =  1.0  x 10~4 
RINHMG8:   Same  as  INHMG8  except 
a  =  0.25  x 10"4 

C.   RESULTS 

This  section  will  present  and  discuss  the  results  from 

the  computer  model  runs  on  the  test  cases  as  defined  in  Section 

IV. B.   Some  cases  were  run  using  both  methods  of  integration 

w.r.t.  u  .   The  two  methods  showed  no  significant  differences, 
o  ^ 

First,  some  trivial  checks  were  made  using  the  various 
homogeneous  cases.   Note,  that  all  test  cases  have  the  same 
line  of  sight  direction  cosine  v  .   For  the  homogeneous  cases 
with  a  single  frequency  component  of  1000  Hz  the  amplitude 
and  phase  terms  of  the  ocean  transfer  function  are  equal  and 
given  by: 


*m   "    <koV"1/2  '       V   =    °  (4-15) 
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Since  those  cases  have  the  same  ocean  transfer  function 
values,  the  difference  in  the  magnitude  of  the  frequency 
response  H(f,i,q)  for  the  various  cases  should  only  depend 
on  the  total  range  given  by 


(x  -x  )2  +  (y  -y  )2  +  (z  -z  )2    (4.16 
total      V   r    o       r  Jo  r    o 


and  the  transmit  far-field  beam  pattern.   Cases  HMGl,  HMG3 

and  HMG4  also  have  the  same  line  of  sight  direction  cosine 

value  u  ,  so  the  influence  of  the  transmit  far-field  beam 
r 

pattern  is  the  same  for  these  cases  and  magnitude  differences 

between  the  frequency  responses  can  only  depend  on  the  total 

range  and  should  behave  in  accordance  with  the  expected 

spherical  spreading  for  a  homogeneous  medium.   Table  IV. 1 

shows  cases  HMGl,  HMG3  and  HMG4  together  with  their  total 

range  and  the  magnitude  of  the  frequency  response  at  the  center 

element  of  the  receive  array,  i.e.,  |H(f,0,0)|.   We  conclude 

that  the  magnitude  indeed  falls  off  as  1/r,  only  the  total 

range  is  of  importance. 

Table  IV. 2  shows  cases  HMGl,  HMG2  and  HMG6  which  differ 

only  in  cross  ranae x  -x   and  line  of  sight  direction  cosine 
J  r    o 

value  u   but  have  the  same  total  range.   In  these  cases  only 
the  far-field  beam  pattern  can  cause  differences  in  the  magni- 
tude of  the  frequency  response.   The  table  shows  that  as  the 
cross  range  increases  the  magnitude  falls  off,  which  is  con- 
sistent with  the  fact  that  when  a  beam  pattern  is  steered 
towards  end-fire  (to  larger  values  of  u,  and  v,)  the 
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directivity  index  decreases  due  to  the  increasing  beamwidth 
of  the  main  lobe. 

Table  IV. 3  shows  the  difference  in  magnitude  of  the  fre- 
quency response  at  the  center  of  the  receive  array  between 
HMG1,  HMG5,  and  HMG7 .   Cases  with  the  same  geometry  but  a 
different  value  for  the  single  frequency  component.   This 
difference  in  magnitude  is  due  to  the  frequency  dependent 
amplitude  factor  inherent  in  the  WKB  approximation,  which  is 
not  exact  for  a  homogeneous  medium. 

Table  IV. 4  compares  the  magnitudes  of  the  frequency 
response  between  the  homogeneous  cases  HMG1,  HMG2  ,  and  HMG5 
and  their  inhomogeneous  counterparts  INHMG1,  INHMG2  and  INHMG5 
We  see  that  the  magnitudes  for  the  inhomogeneous  cases  are 
slightly  lower,  an  expected  phenomena  because  of  the  channel- 
ing effect  in  the  inhomogeneous  medium.   In  a  medium  with  a 
constant  positive  gradient,  the  acoustic  rays  are  bent  upwards 
and  the  sound  intensity  decreases  with  increasing  depth. 
This  causes  a  lower  magnitude  of  the  acoustic  signal  then 
expected  by  spherical  spreading  only. 

We  will  now  present  numerical  and  graphical  output  from 
the  three  dimensional  DFT  beamforming  program.   The  frequency 
transform  for  the  single  frequency  component  cases  simply 
reveals  the  magnitude  of  the  transfer  function  multiplied  by 
the  amplitude  of  the  input  electrical  signal  and  is  not  shown 
here.   Figures  7  through  10  show  the  magnitude  of  the  spatial 
DFT  w.r.t.  x,  i.e.,  Y   (q,r,n)  versus  direction  cosine  u,  and 
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the  magnitude  of  the  spatial  DFT  w.r.t.  y,  i.e.,  Y   (q,m,s) 

versus  direction  cosine  u,  for  the  cases  HMG1,  HMG2 ,  INHMGl 

and  INHMG2 .   These  graphs  allow  us  to  estimate  the  direction 

of  propagation  of  the  incoming  plane  wave  by  finding  the 

coordinates  of  the  maxima  on  the  graphs.   Tables  IV. 5  to 

IV. 8  provide  numerical  values  for  the  above  mentioned  DFT ' s 

w.r.t.  x  and  y.   They  show  values  of  Y   (q,r,n)  and  Y  N(q,m,s) 

around  the  maxima  and  make  it  possible  to  form  more  accurate 

estimates . 

From  the  estimated  values  u   and  v  ,  we  can  calculate 

o      o 

estimates  of  target  elevation  angle  9   and  target  bearing 

angle  (p  .   Table  IV. 9  gives  the  results  for  the  different 

cases.   The  values  of  the  estimators  u   and  v   were  determined 

o       o 

by  performing  a  second  order  interpolation  on  the  data  given 
in  tables  IV. 5  through  IV. 8.   It  is  seen  that  the  estimates 
of  target  angles  for  the  homogeneous  cases  are  correct  and 
equal  the  line  of  sight  angles  from  the  receiver  to  the  trans- 
mitter.  The  estimated  target  angles  for  the  deterministic 
inhomogeneous  cases  deviate  slightly  from  the  line  of  sight 
angles.   This  is  due  to  the  positive  sound  speed  gradient 
which  causes  upward  bending  of  the  acoustic  rays  so  the  direc- 
tion cosine  v  ,  directly  related  to  the  vertical  angle  of 
propagation  of  the  plane  wave,  decreases  as  the  wave  travels 
along.   The  net  effect  for  these  cases  is  relatively  small 
but  noticable. 

The  effect  of  randomness,  characterized  by  a  random  component 
in  the  index  of  refraction  with  a  standard  deviation  of 
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-4 
a  =  1.0  x 10    is  examined  next.   Figure  11  shows  estimates  of 

u   and  v   for  the  case  RINHMG1  and  Fig.  12  shows  estimates  of 
o       o  3 

v   for  cases  RINHMGl  and  RINHMG2 .   As  can  be  seen,  the  esti- 
o 

mate  of  u   is  not  influenced,  which  follows  readily  from  the 
o  J 

fact  that  the  randomness  is  only  present  in  the  y  direction. 

The  estimate  of  v   as  shown  in  Fig.  12  for  the  case  RINHMG2 

o  3 

shows  the  same  random  fluctuations  as  the  one  shown  for  case 
RINHMGl  in  Fig.  11,  since  both  cases  used  the  same  seed  for 
the  random  number  generator,  resulting  in  the  same  random 
number  sequence.   Figure  12  also  shows  RINHMGl  with  a  different 
seed  for  the  random  number  sequence.   Although  both  are  only 
two  possible  samples  they  clearly  show  the  strong  effect  of 
the  randomness  on  the  beamformer  output  and  the  deterioration 
of  the  quality  of  the  estimates  (see  Table  IV. 9). 

Two  three-dimensional  plots  are  shown  in  Fig.  13,  one  for 
the  deterministic  inhomogeneous  case  INHMG2  and  one  for  the 
random  inhomogeneous  case  RINHMG2 .   These  graphs  illustrate 
the  kind  of  output  that  can  be  obtained  from  the  DFT  beam- 
former,  although  it  is  difficult  to  accurately  read  off 
numerical  values. 

Finally,  the  cases  HMG8 ,  INHMG8  and  RINGHM8 ,  which  involve 
several  frequency  components,  are  considered.   The  output 
for  these  cases  is  given  in  numerical  and  graphical  form. 
Table  IV. 10  compares  the  Fourier  series  coefficients  of  the 
output  electrical  signal  with  the  Fourier  series  coefficients 
of  the  transmitted  electrical  signal.   The  widely  varying  values 
are  caused  by: 
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1)  The  beam  pattern  of  the  transmit  array  is  frequency 
dependent.   The  differences  between  the  transmitted  frequency 
components  are  relatively  large  and  cause  large  differences  in 
amplitude  of  the  Fourier  series  coefficients.   The  maximum 
values  are  observed  at  f  =  3000  Hz,  i.e.,  the  frequency  for 
which  beam  steering  was  done. 

2)  The  deterministic  inhomogeneous  medium:   As  already 
shown  in  other  test  cases  the  amplitude  of  the  received  signal 
is  frequency  dependent  due  to  the  WKB  approximation  to  the 
wave  equation  in  an  inhomogeneous  medium.   Equation  (2.10) 
clearly  shows  the  frequency  dependence  of  the  medium  transfer 
function . 

3)  The  random  inhomogeneous  behavior  of  the  medium:   The 
randomness  introduces  additional  deviations.   It  can  be  seen 
from  table  IV. 10  that  the  effect  of  the  randomness  increases 
with  increasing  frequency.   Also,  in  the  case  of  a  random 
medium,  the  estimation  of  the  Fourier  series  coefficients 
becomes  strongly  dependent  on  the  location  of  the  element  in 
the  receive  array.   Note  that  table  IV. 10  only  shows  values 
for  the  receive  element  (m  =  0,  n  =  0) . 

The  received  angular  spectrum  for  the  different  frequency 

components  of  the  cases  INHMG8  and  RINHMG8  is  shown  in  Figs. 

14  through  19.   For  the  deterministic  inhomogeneous  case, 

Figs.  14  through  16  show  a  regular  behavior  of  the  angular 

spectrum  according  to  the  rectangular  amplitude  weighting 

used  for  the  receive  array.   Estimation  of  u   is  exact. 

J  o 
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Estimation  of  v   returns  a  slightly  different  value  from  the 

o  3    J 

line  of  sight  value  for  direction  cosine  v   due  to  the 

inhomogeneous  behavior  of  the  medium.   Figures  14,  15  and 

16  demonstrate  the  increasing  directivity  of  the  receive 

array  with  increasing  frequency,  thus  leading  to  more  accurate 

estimates  for  higher  frequencies. 

Figures  17,  18  and  19  clearly  show  the  strong  influence 

of  the  randomness  on  the  resulting  angular  spectrum.   The 

randomness  creates  strong  "false"  sidelobes  and  causes  an 

offset  of  the  estimated  value  of  direction  cosine  v  .   The 

o 

effects  of  the  randomness  increase  with  increasing  frequency. 

At  f  =  5000  Hz  (Fig.  19)  the  estimate  of  v   is  determined  by 

a  spurious  side  lobe  and  returns  a  totally  wrong  value. 

Table  IV. 11  shows  estimated  values  of  u   and  v   at  the  differ- 

o      o 

ent  frequencies  and  their  associated  estimates  of  target 
elevation  and  bearing  angles,  9   and  ty    ,  respectively. 
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TABLE  IV. 1 

Magnitude  of  Frequency  Response  Comparing  Cases 
Which  Are  Different  Only  in  Total  Range 


Case  Total  Range  |H(f,0,0) 

(in  m) 


HMG1  2000  0.02434 

HMG3  1000  0.04868 

HMG4  3000  0.01622 


TABLE  IV. 2 

Magnitude  of  Frequency  Response  Comparing  Cases 
Which  Are  Different  Only  in  Cross  Range 


Case       Total  Range       Cross  Range       |H(f,0,0 
(in  m)  (in  m) 


HMG1         2000  0.0  0.02434 

HMG2         2000  200.0  0.02417 

HMG6         2000  1224.75         0.01654 
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TABLE  IV. 3 

Magnitude  of  Frequency  Response  Comparing  Cases 
Which  Are  Different  Only  in  Frequency 


Case  Frequency  |H(f,0,0) 

(in  Hz) 


HMG1  1000  0.02434 

HMG5  2000  0.03441 

HMG7  4000  0.04866 


TABLE  IV. 4 

Magnitude  of  Frequency  Response  Comparing 
Homogeneous  and  Inhomogeneous  Cases 


Case  H(f,0,0) 


HMG1  -  INHMG1  0.024  34  -  0.02341 

HMG2  -  INHMG2  0.02417  -  0.02326 

HMG5  -  INHMG5  0.03441  -  0.03310 
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TABLE  IV. 5 

Magnitude  of  Yg^ (q , r , n)  Vs.  Direction 
Cosine  u  for  HMGl,  INHMG1 


Direction 

HMGl 

INHMG1 

Cosine 

u 

0.015 

.0121416 

.0116804 

0.010 

.0121566 

.0116945 

0.005 

.0121655 

.0117030 

0.0 

.0121685 

.0117058 

0.005 

.0121655 

.0117030 

0.010 

.0121566 

.0116945 

0.015 

.0121416 

.0116804 

TABLE  IV. 6 

Magnitude  of  YsN(q/r'n)  Vs  •  Direction 
Cosine  u  for  HMG2 ,  INHMG2 


Direction 

HMG2 

INHMG2 

Cosine 

u 

0.085 

.0120593 

.0116000 

0.090 

.0120744 

.0116147 

0.095 

.0120836 

.0116239 

0.100 

.0120869 

.0116274 

0.105 

.0120842 

.0116254 

0.110 

.0120756 

.0116178 

0.115 

.0120611 

.0116046 
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TABLE  IV. 7 

Magnitude  of  YsN(q/m/s)  Vs •  Direction 
Cosine  v  for  HMGl ,  HMG2 


Direction 

HMG1 

HMG2 

Cosine 

V 

0.484 

.0121367 

.0120534 

0.489 

.0121531 

.0120696 

0.494 

.0121636 

.0120800 

0.499 

.0121681 

.0120844 

0.504 

.0121666 

.0120828 

0.509 

.0121591 

.0120754 

0.514 

.0121457 

.0120620 

TABLE  IV. 8 

Magnitude  of  YSN(q,m,s)  Vs.  Direction 
Cosine  v  for  INHMGl,  INHMG2 


Direction 

INHMGl 

INHMG2 

Cosine 

V 

0.475 

.0116731 

.0115947 

0.480 

.0116893 

.0116111 

0.485 

.0117001 

.0116215 

0.490 

.0117051 

.0116264 

0.495 

.0117045 

.0116258 

0.500 

.0116983 

.0116196 

0.505 

.0116864 

.0126078 
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TABLE  IV.  9 


Estimated  Values  u  ,  v  ,  0   and  i> 

o    o    o      o 


Case 

u 
o 

V 

o 

(in  degrees) 

HMG1 

0.000 

0.500 

+30.00 

HMG2 

0.100 

0.500 

+30.65 

INHMG1 

0.000 

0.492 

+29.47 

INHMG2 

0.101 

0.492 

+30.15 

RINHMG2 

0.101 

0.523 

+32.19 

RINHMG1 

0.000 

0.523 

+31.53 

RINHMG1 

0.000 

0.508 

+30.53 

(in  degrees) 

+  90  .00 

+78.69 

+90.00 

+78.40 

+79.07   (seed  #1) 

+90.00   (seed  #1) 

+90.00   (seed  #2) 


1)   The  line  of  sight  angles  from  receive  to  transmit 
array  (or  target  angles)  are: 

HMGl,  INHMG1,  RINHMGl :   9   =  +30.00°,  ip      =    +90.00° 


HMG2,  INHMG2,  RINHMG2  :   6   =  +30.65°,  ip   =  +7  8.69 


2)   Note  that  equation  (4.14)  which  is  used  to  calculate 
estimates  of  target  bearing  angle  from  estimates  of 
uQ  and  v0  returns  two  results.   This  is  due  to  the 
ambiguity  of  a  planar  array.   For  the  above  table  the 
appropriate  values  were  selected. 
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TABLE  IV.  10 

Estimation  of  Fourier  Series  Coefficients 
YSN(q,0,0)   for  HMG8,  INHMG8  and  RINHMG8 


Frequency 

HMG8 

INHMG8 

RINHMG8 

(in 

Hz) 

1000 

0.0010987 

0.0009705 

0.0010341 

2000 

0.0098126 

0.0098901 

0.0102369 

3000 

0.0209366 

0.0201906 

0.0177316 

4000 

0.0138747 

0.0122052 

0.0087726 

5000 

0.0024539 

0.0033176 

0  .0051316 

Note : 

The 

Fourier  series 

coefficients 

of  the  trans- 

mitted 

electrical  signal  are  c. 

=  ak  =  1; 

ri     "~"   J_/<--  *«••/_> 


TABLE  IV. 11 


Estimated 

Values 

of  u  , 
o 

v  ,  9   and 
o    o 

^o 

/S 

/s 

/N 

/\ 

Case 

Frequency 

u 
o 

V 

o 

eo 

^o 

(in  Hz) 
1000 

(in  degrees) 
+  31 

(in  degrees) 

HMG8 

0.10 

0.50 

+  79 

ii 

2000 

0.10 

0.50 

+  31 

+  79 

it 

3000 

0.10 

0.50 

+  31 

+  79 

ii 

4000 

0.10 

0.50 

+  31 

+  79 

ii 

5000 

0.10 

0.50 

+  31 

+  79 

INHMG8 

1000 

0.10 

0.49 

+  30 

+  78 

ii 

2000 

0.10 

0.49 

+  30 

+  78 

ii 

3000 

0.10 

0.49 

+  30 

+  78 

ii 

4000 

0.10 

0.49 

+  30 

+  78 

M 

5000 

0.10 

0.49 

+  30 

+  78 

RINHMG8 

1000 

0.10 

0.63 

+  40 

+  81 

ii 

2000 

0.10 

0.49 

+  30 

+  78 

ii 

3000 

0.10 

0.53 

+  33 

+  79 

ii 

4000 

0.10 

0.52 

+  32 

+  79 

ii 

5000 

0  .10 

-0.43 

-26 

-77 
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Fig.  7    Estimation  of  Direction  of  Propagation 
for  Case  HMGl 
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Fig.  8    Estimation  of  Direction  of 
Propagation  for  Case  HMG2 
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Fig.  10    Estimation  of  Direction  of 
Propagation  for  Case  INHMG2 
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Fig.  11   Estimation  of  Direction  of  Propagation  for 
Case  RINHMG1 
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Fig.  12   Estimation  of  v  for  Cases  RINHMG2  and 
RINHMG1  (Different  Seeds) 
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Fig.  15   Angular  Spectrum  for  Case 
INHMG8  at  f  =  3000  Hz 
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Fig.  19   Angular  Spectrum  for  Case 
RINHMG8  at  f  =  5000  Hz 
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V.   CONCLUSIONS  AND  RECOMMENDATIONS 

Computer  simulation  of  a  mathematical  model  of  wave  propa- 
gation through  a  random,  space-variant,  time-invariant,  ocean 
medium  between  two  planar  arrays  has  been  accomplished.   The 
computer  simulation  incorporates  an  index  of  refraction  which 
is  only  a  function  of  depth.   The  deterministic  component  of 
the  index  of  refraction  represents  a  sound  speed  profile  with 
constant  gradient.   The  random  component  is  assumed  to  have 
statistical  properties  independent  of  depth.   The  effects  of 
the  randomness  on  the  electrical  output  signals  at  the  differ- 
ent receiver  elements  are  assumed  to  be  uncorrelated. 

The  computer  program  was  written  in  FORTRAN.   It  implements 
both  analytical  and  numerical  integration  techniques.   The 
results  from  both  techniques  are  in  close  agreement. 

A  number  of  test  cases  have  been  run.   Interpretation  of 
the  data  for  simple  cases  show  results  which  are  consistent 
with  expectations  from  theory.   The  model  correctly  demonstrated 
the  influences  of  the  far  field  beam  pattern  of  the  transmit 
array,  the  source-receiver  geometry,  and  the  deterministic 
inhomogeneous  behavior  of  the  medium.   Also,  a  random  inhomo- 
geneous  medium  was  shown  to  affect  the  wave  propagation  and 
distort  the  frequency  and  angular  spectra,  leading  to  errors 
in  the  estimation  of  the  Fourier  series  coefficients,  target 
bearing  angle,  and  target  elevation  angle.   A  final  test  case 
RINHMG8,  with  a  transmit  signal  containing  several  frequency 
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components  demonstrated  a  complex  dependence  on  the  above 
mentioned  factors,  making  interpretation  of  the  results  diffi- 
cult.  The  effect  of  the  randomness  was  very  pronounced  in 
this  case. 

The  computer  program  could  be  further  improved  by  imple- 
menting a  more  sophisticated  method  to  simulate  the  randomness 
of  the  medium.   The  assumption  that  the  random  phase  terms  at 
the  different  receiver  elements  are  uncorrelated  breaks  down 
for  small  interelement  spacings.   The  computer  program  can 
also  be  easily  modified  to  accommodate  a  more  complex  sound 
speed  profile  in  which  the  gradient  is  a  function  of  depth, 
for  example. 

Future  application  of  the  model  could  be  devoted  to 
investigating  the  propagation  of  more  complex  transmit  signals 
One  could  also  investigate  the  effects  of  the  randomness  of 
the  ocean  medium  on  signal  processing  using  different  sizes 
of  receive  arrays,  i.e.,  increasing  the  number  of  elements. 
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